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I. Introduction 


The objective of the study reported herein is to infer the exterior 
gravitational field of the earth and the structure of the geoid from GE(^S-C 
metric data. Of the many measurement types available, the scope of our 
effort is limited to two: satellite-to- satellite range and range rate tracking 

(SST) and radar altimetry. Initially, our intent was to couple together the 
dynamical effects of the field perturbations on doppler or range rate with 
the geopotential- induced geoid displacements sensed by altimetry to obtain a 
single consistent model. However, it soon became apparent that the coupling 
between the two data types was too weak to justify the complexity involved 
in the short arc method of analysis employed. Because of the weak coupling, 
the analysis effort separated naturally into two parts: inference of the long 
wavelength disturbances from SST and the fine structure geoid from altimetry. 

The division was also fortuitous from a chronological standpoint. The 
SST tracking which was actually used occurred prior to June 28, 1975 and all 
of it became available by mid-1976. Many passes of altimetry were also avail 
able but the coverage over the test area was not sufficiently dense until the 
latter part of 1976. By that time the SST work was well underway. 

Section II of the report concerns various phases of the SST analysis. 

A direct point- by-point estimate of gravity disturbances by means of a recur- 
sive filter with backward smoothing was initially attempted but had to be 
forsaken due to poor convergence (attributed to inadequate signal-to-noise 
ratio). The adopted representation consists of a more or less uniform grid 
of discrete masses at a depth of approximately 400 km from the earth's 
surface. The layer is superimposed on a spherical harmonics model, in 
this case GEM- 6 (Ref. 1). 

The second part of the report describes the procedure for smoothing 
the altimetry and inferring the fine-structured gravity field over the Atlantic 
test area. In place of the usual surface anomalies, we again chose to repre- 
sent the local disturbances by means of a density layer. The altimeter height 
biases were first estimated by a least squares adjustment at orbital cros sover 
points. The bias-adjusted data base contained 74 arcs with 6050 observations 
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and 1150 crossovers. After taking out the bias, long wavelength contribu- 
tions from GEM-6 as well as a calibration correction from Reference 2 
were subtracted. The residual heights were then represented by a mass 
distribution beneath the earth's surface. 


II, Analysis of Satellite-To-Satellite (SST) Tracking 


The geometry of SST tracking is illustrated in Fig. 1, ATS-7 at 
synchronous altitude serves as a link between ground station and GEOS-C 
for range and doppler measurements. Two basically different approaches 
of analysis can be employed depending on the length of the data span. Some 
general comparative features are: 

Long Arcs Short Arcs 

Methods Methods 

a) Length of orbit/arc Many revolutions Less than one revolution 

(-^100) /arc per arc 

b) Number of arcs Few - 3-4 over Many - 50-60 over 2 months 

several months 

c) Size of initial state Small - 20-30 Large - 350 

vector 

d) Aliasing Maybe significant A daptable to uneven data 

when data distribution distribution, some absorp- 

is uneven tion of initial perturbations 

by state vectors 

e) Sensitivity More sensitive, Sacrifices sensitivity to 

particularly to long long wavelengths 

period perturbation 

f) Orbital dynamics Favors dynamical Numerical integration 

theory for long term applicable 

integration 

g) Dynamical model Orbital theory avail- Harmonics, discrete 

able for harmonics, anomalies, and all modeled 

planets, solar pressure forces 

h) Linearity May be nonlinear Very linear, usually con- 

unless one starts with verges in a single iteration 

good state vector with reasonably good initial 

states 

i) Singularity Nonsingular in states, May be near singular in both 

may be near singular states and gravitational param- 
in harmonics eters. Some a priori informa- 

tion stabilization usually re- 
quired depending on coverage 
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Figure 1. Geometry of SST Tracking 
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j) Computer time Not evaluated 

requirements/ 
iteration 


k) Accuracy Adequate 


Acceptable - 3200 sec/ 

iteration on CDC 7600 for 
71 6 parameters, 4310 
observations; large part 
of time taken up in computing 
observations and observation 
partials 

Adequate down to less than 
. 1 mm/ sec 
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We chose the short arc approach for reasons of d, f, g, h, j above 
and also because of our past experience with the method in analyzing lunar 
and Mars orbiter data. Our inventory of software required relatively minor 
revisions. To summarize, the high degree of linearity, reduction of process 
noise (aliasing), numerical accuracy, simplicity, computer time economy, 
and software availability were the factors in our choice. 


Having decided on a short arc technique, there are a number of consider- 
ations such as the filtering algorthm, data types to analyze, and the method of 
representation. Initially our intent was to combine SST with ground tracking, 
i. e. , doppler, C-band, laser, together with altimetry to obtain a single con- 
sistent solution. Each element of the disturbing potential, T, would produce 
an acceleration VT on the motion of GEOS-C sensed by the various tracking 
systems. It would also cause a height anomaly, AN = T/Y , measurable by 


altimetry and a gravity anomaly, Ag 


derivable from surface 


_ ^ 

dr T 

data. The plan is ideal in principle and might be amenable to analysis utilizing 
a long arc approach. For short arcs, the combination of the data types is 
unnecessarily complicated. Since GEOS is at an altitude of 850 km, one would 
not expect its orbit to be sensitive (over time spans of a half revolution) to 
disturbances with wavelengths, much less than twice the altitude. Since the 
orbit perturbations are sensing features with \ on the order 2000 km, pertur- 
bations arising from fine structure variations of few tens kilometers as 
measured by the altimeter would be entirely lost in the noise. There was 
no need to add the altimetry to the SST data in a short arc analysis. Decoupling 
simplified the problem considerably. 


There is next the question of a filter or smoothing algorithm for the SST* 
Some possible choices are: 1) arc by arc polynomial spline fits of doppler data 

residuals which are based on a low degree harmonics model followed by analytical 
differentiation to obtain point accelerations, 2) point by point estimate of un- 
modeled accelerations with a recursive filter--also one arc at a time, 3) directly 
fitting to all selected arcs by means of a set of discrete parameters (may be 
harmonics coefficients, density elements, surface anomalies, or geoid heights). 


•■'X 

Y is normal gravity and r is the radial coordinate. 
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The Jet Propulsion Laboratory has applied the first extensively to the tracking 
of lunar, orbiters, Apollo and GEOS .[l, 2, 3]. Only limited application of the 
second has been made [4, 5]. The last is the one used on the Moon and Mars 
[6, 7 ]. The principal differences between the first two and the third is that, in the 
former, accelerations are derived directly from the data and the sources can 
be inferred by means of a static (nondynamical) fit at a later time. The 
advantage of the third method over the first is that aliasing, where a part of the 
gravitational signal residuals may be absorbed by the states in the initial fit 
with a low degree model, is reduced by estimating state and field parameters 
simultaneously. In principle, method 2 is the most attractive because it com- 
bines simultaneous state and unmodeled-acceleration estimation with the 
possibility of static fit. The latter feature is so important that we decided to 
examine the feasibility of filtering in some detail by means of a simulation study 
for possible application to the analysis of the SST tracking. 

2.1 Recursive Filtering Simulations 

The steps in the filter simulations were as follows: 

a) Generate range and range rate observations from ATS and GEOS 
using a dynamical model containing small gravitational distur- 
bances (on the order of 1-10 milligals at GEOS altitude) super- 
imposed on a low degree field. Zero means Gaussian noise 
with standard deviation of 1 meter and 1 mm/sec, respectively, 
were added to the generated data. 

b) A forward and backward (or smoothing [s]) recursive filter with 
plant noise compensation was applied to the generated data to 
estimate the unmodeled accelerations. The filter dynamics 
contained the nominal harmonic field plus the filtered acceler- 
ations. At each time point, t, the forward and backward filter 
outputs were combined to obtain a minimum variance estimate 
including the effects of both data before and after t^. 
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Differential equations for the filter are: 


x = v; v=a4-u; u= - pu + 0 

where x, v, a, and u are 3-vectors of position, velocity, modeled acceleration, 
and unmodeled acceleration, respectively. 0 is a zero mean Gaussian process 
with given covariance, p is a reciprocal time constant. Parameters varied 
in the simulation included a priori variances on the states, p \ plant noise 
level. A wide range on the parameters were tried: 

Parameter 
cr(v) 

<r(Uo) 

p-1 

<r(0) 

Initially, the unmodeled acceleration was limited to about 2 mgals maximum 

amplitude. The filter results were disappointing. When(r(0) was very small, 

-2 

say 10 mgals, (approximate random part of the field over a 10 second span) 
a posteriori <t(u) was on the order of . 05 - . 1 mgals for values of P ^ ^ 500 
seconds, (p was held fixed). Increasing cr(0) to .2 mgals increased <y(u) to a 
range of . 5 - 3 rngals. The latter is certain to be too large. In either case, 
the filter output deviated from the correct answer by an unacceptably large 
amount--from a few tenths to as much as 10 mgals. The results were not 
altogether unreasonable since two observations over a 10 second interval 

2 

could, at best, measure acceleration along the line-of-sight to say . 14 mm/ sec 
— 14 mgals. With such a low signal to noise ratio, a lengthy smoothing time 
(or large value of p ^) was needed and the random component cannot be assumed 
correlated over times much greater than, say, 100 seconds. An experiment 
in which the perturbation level was increased ten fold (to a maximum of 20 mgals) 
was also made. A typical output is shown in Fig. 2a which should be compared 


Range 

5 - 1000 meters 
. 005 - 1 m/sec 
1-10 mgals 

loo - 1000 seconds 

_2 

10 - 3 mgals over 10-40 second interval 
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with the true values shown in Fig. 2b. The figures have some resemblence 
to each other, but are clearly too far apart for the level of accuracy being 
sought. The recursive method would undoubtedly be more appropriate for the 
Moon or Mars where anomalies on the order of several hundred milligals are 
present. In the case of the earth, a, measurement noise level of 1 mm/ sec 
appears too high for detecting the small signals at GEOS altitudes by filtering. 

2. 2 Inference Procedure 

Having forsaken the filtering approach, we reverted to the method of 
representation by discrete eleinents, in this case point masses. The object 
is to infer a mass or density layer by fitting to the doppler observations. In 
the case of SST tracking under consideration, an outline of the steps in the 
inference procedure is as follows: 

a) Select the data base for ATS and GEOS. Edit and presmooth 
data, separate into individual passes of about one hour duration. 

b) Determine accurate ATS ephemeris from data base and update 
to GEOS tracking periods. 

c) Obtain preliminary GEOS states at beginning of each pass. 

Adjust states to obtain best fit of SST tracking using model 
which includes all significant forces except residual gravity 
disturbances. Each arc is adjusted separately. 

d) Establish grid of mass point coordinates. Generate partial 
derivatives of doppler observations with respect to mass 
magnitudes and the six initial states. Starting with those from 
step c, the partials are generated one arc at a time. 

e) Form nox-mal equations and solve to find improved states and 
masses under constraints. 

f) Iterate to convergence. 

2. 3 Data Preparation 

Much of the prelimiirary work was related to decoding, presmoothing, 
editing and reformatting of the data tape. A program was written to perform 
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all four functions and punch cards in a format acceptable by the Aerospace 
orbit program as output. In addition to station code, time, and smoothed 
measurement, each card also contained frequency synthesizer mode codes 
to enable calculation of appropriate scaling constants for the doppler count. 
Approximately 7500 smoothed observations were selected, 3000 ATS and 
4500 GEOS, representing 12000 and 18000 uncompressed 10 second samples 
for ATS and GEOS, respectively. 

Pre -Smoothing and Editing 

Computation of the observable, i. e. , the doppler count and its partial 
derivatives, consumed a substantial part of the time for the regression because 
of the tedious iterative algorithm involving five range legs at the beginning 
and end of each count interval. To conserve time, each frame of four 10 second 
points was compressed into a single smoothed point by means of a quadratic 
polynomial. A cubic could also be used but the cubic yields exactly the same 
value at the midpoint as a quadratic with the difference that the residuals are 
zero. The quadratic was used since it was simpler and because the residuals 
from the second degree fit provided a means for the initial editing of outliers 
(of which there were a significant number, particularly near the extremities 
of a pass). In principle, the smoothed values should have a smaller noise 
variance but the reduction is quite small- -on the order of 20 percent for the 
standard deviation. The main purpose is that the 40 second sampling rate 
reduced computing time by a factor of four with hardly any sacrifice in 
Sensitivity. 

Refraction induced residual tails were very much in evidence for all the 
long passes where the ray path traversed the ionosphere. Since modeling of 
the refraction profile was not one of our objectives, the first and last few 
minutes of data in each pass were dropped. The earth central angle, 0, 
between ATS and GEOS served as a quantitative rejection criterion. Obser- 
vations for which 0 exceeded 95. 5 degrees were deleted. The cutoff value 
corresponds to a minimum altitude of 600 km for the ray path. 

A similar data compression scheme was applied to ATS. The loss of 
sensitivity is even less in the case of ATS. 
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Data Base ^ GEOS 


The original data tapes contained about 176 passes taken over the period 
from 16 April to 26 June 1975. Since the GEOS altitude was ~ 850 km, we 
wanted to select passes with an equatorial longitude spacing on the order of 
the altitude or slightly less than 8 degrees. Forty five passes would cover 
the entire circumference and 90 would suffice to include both ascending and 
descending geometries. The actual passes fell far short of the desired 
coverage. Much of the Atlantic had redundant passes while the Pacific was 
not included at all (such an uneven distribution tends to favor a discrete or 
mass layer representation over a spherical harmonics expansion). In the 
process of selecting the passes to be used in our solution, we examined 
nearly every pass available as to quality, geometry, and coverage. A total 
of 59 passes were included iii the final data base. Figure 3 shows the geographic 
distribution of the pass while Table 1 summarizes some pertinent statistics. 

The average number of points (nominally 40 seconds apart) per pass is 73 
and the average duration about 50 minutes. 

The pass geometry may be graded by the minimum value of 6, Passes 

with - 50 degrees comprise about 75 percent (44 passes). The remainder 

have 0 . between 50 and 65 degrees. It should be mentioned that 0 ranges 

mm a to 

from about 95. 5 degrees to 0 . on all the passes, i. e. , 0 is at or near 0 

^ min ^ min 

(the vicinity where the doppler data has the maximum sensitivity to anomalies) 
only over a fraction of the pass (see next section). 

The coverage of the earth by good geometry passes is far from complete, 
even for the regions from -lOO to +100 longitude --but we are unable to fill in 
any more gaps from the data tape currently available. Many of the candidate 
passes in early May either were too noisy or had dropouts during the crucial 
segments when the sensitivity is near the peak. 

2. 4 Software Features and Modifications 

The principal analysis tool is a, modified version of the Aerospace orbit 
program, TRACE. The program numerically integrates the dynamical 
equations, including effects of sun, moon, geopotential, and radiation pressure 
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6 

c c 

6 

11 

7 

3 

37 

71. 3C 

A 

99 

71 

1037 

6 

1 Z 


4? 

8 

40 

17 

46 • C C 

A 

c c 

7C 

1057 

6 

cl 

17 

11 

17 

51 

50 

92. 2 C 

ca 

c c 

e c 

0 

1063 

6 

24 

4 

8 

5 

0 

60 

107. 6C 

AD 

99 

73 

1064 

6 

24 

'C 

4» 

6 

32 

43 

32. 2C 

A 

9-9 

€4 

1078 

6 

25 

c 

31 

6 

19 

46 

87.70 

A 

CC 

64 

10 85 

6 

2 S 

16 

35 

17 

24 

55 

103. 1C 

ca 

Is 

6G 

1091 

6 

2 6 

3 

42 

4 

33 

62 

118. 6C 

ao 

CC 

69 

1092 

6 

2 6 

C 
^ . 

14 

6 

C 

49 

93. 2C 

a 

CC 
0 0 

65 


*YEAR = 1975 , MN = month, DA = day, HR = hour, THETMIN = Minimum 
earth central angle from ATS to GE(JS, NODE = east longitude in degrees; 
AS, A = ascending, D = descending; ST, 58 = Rosman, 99 - Madrid, NPT = 
number of 40 second samples. 
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(for ATS only) in an inertial coordinate system, as well as the variational 
equations for parameters being estimated. Both trajectory and variational 
equation solutions are interpolated at observation times which are intermediate 
to integration step times. A subroutine then computes the value of the obser- 
vations and partials from the interpolated states for the differential correction 
routine. Some relevant properties of the program for this application are: 

a) Coordinate system - Integrates in inertial coordinates; equator, 
equinox at midnight of date, outputs in equator and equinox of date. 

b) Time - Integration time is uniform, output time is UTC, linear 
equations relate integration time to UTC, UTl, and ephemeris 
time. 

c) Geopotential - Includes dynamical and variational equations for 
spherical harmonics as well as discrete masses. 

Some pertinent constants for the analysis are: 

Earth model for GEOS - Full GEM-6 with 174 pairs of terms. Ref [9] 

ATS - GEM-6 truncated at degree 6 

Earth radius - 6378145 meters 

Flattening - (298. 255)”^ 

GM - 3.986008 X 10^^ m^/sec^ 

Velocity of light - 2. 997924562 x 10® m/sec NBS(1972) 

Station coordinates - Best available locations from Mr. J. Marsh 

of GSFC 

Major modifications for the present analysis task included: 

a) Range and range rate equations for satellite-to-satellite tracking 

b) Range trilateration equation for ATS 

c) Variational equations for SST with respect to geopotential 
parameters which had to be expanded to accommodate 362 
point masses plus 6 initial states 
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d) Frequency synthesizer equations to calculate constants for 
scaling the range differences into doppler count. 

A considerable part of the TRACE work involves expanding of sub- 
routines to accommodate the integration of some 362 sets of variational 
equations for the mass grid. The range and range rate algorithms, frequency 
synthesizer equations are given in informal notes from NASA Ref. [10]. 
Ambiguity resolution of both the doppler count and range rollover are done 
in TRACE with the aid of an initial orbit. 

2. 5 ATS Orbit Determination 

Even though the preliminary analysis indicated that the ATS status could 
be observed in idealized geometry cases from a single pass of SST data, the 
fact that ATS state errors of several kilometers and lO's of mm/sec had little 
effect on the gravity estimates prompted us to fix the ATS states in the gravity 
solution thus eliminating a potential source of singularity--particularly in 
cases with degraded geometry. It was inconvenient to use the NASA furnished 
ATS initial states due to probable differences in coordinate system, time, 
and dynamical model. Hence, it was decided to determine the ATS orbits 
from the observations supplied. The results would also serve as a check on 
our algorithms and software. 

The tracking data on ATS from 4/27 to 6/28/75 were divided into seven 

arcs whose length varied from about 3 to 10 days. Care was taken against 

having a maneuver within any segment. We estimated six state parameters, 

a radiation pressure constant, Y , and a number of range bias parameters 

in each arc to fit range and range rate. In contrast to the gravity inference, 

the range measurements contribute much more towards establishing the 

-7 2 

ephemeris than the doppler. An average value of Y = . 727 x 10 kg/m was 
obtained over the two month period. The number of bias parameters depended 
on the data included in the span. A parameter was introduced whenever one 
was indicated either from the tracking notes sent with the ATS tape or from 
the residual patterns. Some statistics on the ATS ephemeris determination are 
given in Table 2. 


16 





Table 2. 

ATS Data 

Fit Statistics 


Start T 

ime 

Stop T 

ime 

Stations ' 

RMS Residuals 

Month 

Day 

Month 

Day 


Range - m 

Doppler - hj 

4 

27 

5 

4 

R, S, Me 

2.9 (245)*'^" 

.010 ( 35) 

5 

11 

5 

20 

R,Me 

11. 2 (724) 

. 0048 (384) 

5 

22 

5 

30 

R, Me, Md 

14. 8 (618) 

. 0086 (324) 

5 

31 

6 

8 

R,Md 

19. 1 (381) 

. 023 (239) 

6 

1 1 

6 

20 

A, Md 

1 1. 1 (373) 

.0059 (245) 

6 

21 

6 

24 

A, Md 

8. 5 (285) 

. 020 (139) 

6 

25 

6 

28 

A,Md 

6. 8 (299) 

. 0059 ( 46) 


= Rosman, Me = Mojave, A = Ahmedabad, Md = Madrid 
Number in parenthesis = number of smoothed observations 

Systematic trends were evident and the fits could be improved by adding 
more bias parameters. In one experiment on the first arc, we achieved an 
RMS of . 8 meters over 1032 range measurements by introducing 34 bias 
parameters. Assuming a measurement standard error of 1 m on each obser- 
vation, the formal statistics indicate 1 sigma uncertainties of; 

<7 (Range bias) 15 meters 
(^(Position) 150 meters 

<7 (Velocity) 1 cm/ sec 

The converged state vector difference between the orbit with the many bias 
parameters (and small RMS residual) and the one with just a few biases and 
poorer fit was about 1 km in position and 5 cm/ sec in velocity. It was felt 
that the biases remaining did not alter the ATS ephemeris enough to affect the 
short arc solutions. Hence, the later arcs had fewer bias parameters in the 
estimations and poorer RMS fit. 
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After the solutions converged, the ATS ephemeris had to be updated 
to epochs just prior to the start of each GEOS arc. They were subsequently 
held fixed for analysis of the SST data. Cornparisions of our ephemeris with 
the ones given in the NASA tracking notes indicate differences of a few kilometers 
and several cm/ sec which should be adequate. Preliminary covariance analysis 
indicates that the ATS errors resulting fronr estimation along with the GEOS 
states in a short arc mode using the SST tracking would be on the same order. 

Z.6 GEOS States 

Considerable effort was expended in prelmiinary experiments with GEOS 
orbit determination. Most of the trial runs relate to combining data from 
various sensors, i.e. , C-band, laser, and SST. The objective was to see 
how much ground tracking contributes in improving observability of orbit and 
gravity field parameters. Since our plan included analysis of altimetry data at 
the time, the accuracy of GEOS altitude was also of particular interest. A 
summary of an example pass on 27 April 1975 is given in Table 3. This pass 
contained some simultaneous C-band, laser, and SST tracking. The latter 
extended from 29670 to 33010 sec GMT and was classed as good pass 
geometry (0 . =^30 degrees, cf Sec. 2. 3). In this run we included C-band 

tracking from a prior rev as well as the one coincident with the SST. Both 
C-band and laser were of excellent quality. Typical laser residuals are shown 
in Fig. 4, The RMS residuals were as low as . 15 m. Formal statistics from 
various data combinations listed in Table 3 revealed some pertinent and more 
or less expected results: 

a) Ground tracking is better than SST doppler in determining GEOS 
ephemeris even if one assumed perfect ATS ephemeris along 
with the latter data. 

b) If one combined ground tracMng with SST and estimated ATS 
and GEOS states, then the GEOS state errors are about the 
same as the ones in case 3, i. e. , ground data alone, but the 
ATS ephemeris is no better than about 1600 meters. Hence, 
little is gained by trying to improve the ATS orbit with ground 
and SST tracking on GEOS. 
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Table 3. Formal Statistics for Sample Case of GEOS Orbit Determination 


Initial GEOS State 

Nominal 

Standard Deviation (a) of GEOS States 

ATS States 

i 

Epoch = 1975, April 27, 7^ 

States 

Case 1 

Case 2 

Case 3 

Case 4 

Case 5 

Case 2 

Case 5 

Right ascension (Ol) - deg 

281. 092 

. 15E-5 

. 44E-5 

. 44E-5 

.38E-4 

. 0052 

. 0013 

. 0039 

Declination (i) - deg 

10.751 

. 20E-5 

. 65E-5 

. 70E-5 

. llE-3 

. 0033 

. 0022 

. 024 

Velocity angle from vertical (^) - deg 

90. 023 

. 23E-5 

. 12E-4 

. 15E-4 

. 30E-5 

. 00010 

. 0030 

. 020 

Inertial azimuth (A) - deg 

334. 541 

. 28E-5 

.90E-5 

. 90E-5 

. llE-3 

. 0028 

.0017 

. 019 

Radius (r) - meters 

7219705. 

. 12 

00 

o 

. 96 

.64 

69 

1580 

9040 

Velocity (v) - m/sec 

7433. 811 

. 120E-3 

. 82E-3 

. 99E-3 

. 59E-3 

.050 

. 039 

. 16 

Case 

Data in Fit 


Regression Parameters 

Rank 



1 

CB, L, SST 


G + 8 biases 


6 



2 

CB, L, SST 


O 

4- 

00 

biases 


12 



w. .. , 3 

CB, L. 


G -f 8 biases 


6 



kO 

4 

SST 


G 



6 



5 , 

SST 


G + A 



9 




Notes; 

1. CB = 2 C'-band trackers, L = 2 laser stations, SST = doppler only (too few range points which also appeared 
biased). 

2. G = 6-vector of GEOS states, A = 6-vector of ATS, 8 biases are range, azimuth, elevation for each C-band 
and range for each laser. 

3. Pass was on 1975, April 28, 8^ 14^ to 9*^ 10^, rev 245/246. 



c) SST data by itself should not be used to simultaneously estimate 
ATS and GEOS ephemeris. Uncertainties in ATS of almost 2 
km are incurred even with a reduced rank solution, in this case 
rank = 9. 

While these runs showed the importance of ground tracking in establishing 
the ephemeris, it was clear that the levels of accuracy obtainable with the SST 
data alone was adequate for sensing the gravity disturbance field. Also, 
previous covariance analysis indicated that little is to be gained from the 
Im range tracking in sensing the small gravity disturbances. Very accurate 
doppler tracking from, the gound would be of use but a search of the data catalog 
revealed very limited simultaneous tracks with SST. Furthermore, the tracks 
are limited in time duration and coverage. In view of the small return for 
labor expended, we decided to use only SST tracking for the gravity solution 
with fixed ATS states derived from ground based range tracking. This procedure 
simplifies the (problem by eliminating some 59 x 6 = 354 ATS state parameters 
from the linear system. 

The GEOS states we started with were obtained from converged orbit 
determination fits to the SST doppler using TRACE with sun, moon, and full 
GEM-6 geopotential. Each arc of about 50 minutes determines a separate set 
of states. Starting values came from the NASA tracking notes extended to the 
initial epoch of each arc, also using the GEM-6 model. These initial runs 
served to edit outliers and the pass tails affected by refraction. Passes with 
biases, discontinuities, poor or redundant geometry were thus eliminated. 

The small remaining residuals from these runs indicated a low signal to noise 
ratio with the GEM-6 model as base harmonics. Reruns on several passes 
using SAG 73 were made. The signatures were higher but still quite low. 

The excellent fit of the data with GEM-6 is a mixed blessing. The small 
magnitude of gravity error means that there will be very little aliasing between 
initial state and gravity model error. It also means that the corrections to 
the model and state will be highly linear because of small amplitudes. On the 
minus side, it means that we will be hardpressed to improve on the model 
unless the observation noise level is low enough to reveal some definitive 

’ ( 
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signature. All of our preliminary orbit determination runs show that the 
best we can hope for is about a 2 to 1 reduction in RMS and the final runs 
show that it is considerably less than 2 to 1 . 

Most of the runs converged in 2 to 3 iterations from the initial states 
with RMS doppler residuals of about .01 Hz. Typical statistics are given in 
case 4 of Table 3, 

2. 7 Mass Distribution Location 

In using a mass layer representation, choices must be made on the 
density of grid points, geographic location, and the extent of the layer if non- 
global. We wanted to extend the distribution such that the distances from 
GEOS to the most remote mass is at least 5 times the distance to the closest 
mass (which is in the vicinity of the subsatellite point) so that the acceleration 
magnitude for the former is reduced by a factor of about 25. At an altitude 
of 850 km the lateral distance should be about 40 deg. From the coverage 
shown in Fig. 1, with the 40 deg extension on each side, there is little un- 
covered area remaining. On this basis, we decided that the layer may as 
well extend over the entire earth when all the arcs are analyzed simultaneously. 
A complete layer allows the easy imposition of certain linear constraints to • 
the density solution. 

Our rule-of-thumb for grid point separation is that it should be about 
equal to the vertical distance to GEOS. If the masses are located on the surface, 
the separation should be about 850 km. The initial thought was to use 10° x 10° 
equal area squares which number 4l0 with a spacing of 1110 km. However, 
the grid has some asymmetries which could cause distortions when the low- 
degree harmonic constraints are applied (see Sec. 2, 8). The actual grid 
chosen is a subdivision, of an icosohedron with a total of 362 points. Each of 
the 20 faces has an identical grid point pattern. The average distance between 
points is ~1270 km. By choosing the radius of the layer to be 5950 km (depth 
425 km), the mean altitude to GEOS is also about 1270 km. Table 4 which 
contains the mass values from the final solution also lists all of the grid point 
latitudes and longitudes. The 10 reference vertices are located at the north 
and south pole; 26. 565°N, n72°E , 26. 565°S, (n72 + 36)°E; where n = 0 to 4. 
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Fig. 4b. Typical Jjaser Range Residuals 
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Normal Equations, Constraints 


In this section we use the notations: 

O - n vector of observations; n = 4310 

X - m vector of parameters to be estimated; m = 716 

W - Matrix of observation weights, generally assigned to 
be reciprocal variance of each observation 

90 

A - nx m matrix of partials, A , A' = A transpose 

9x 

b - n vector of residuals (observed minus calculated) 

In order to minimize the quadratic sum S = b'Wb, we solve the linear system 

N = (A'WA +r"^) X = A'Wb (1) 


where r is the covariance of the initial estimate, Xq. The new estimate is 

+ X and the covariance matrix is (A'WA +F ) . The matrix N is 

accumulated and stored as the submatrices A^ for the i^^ arc are generated 
individually. The vector x consists of 354 initial states (6 each for 59 arcs) 
and 362 mass parameters. 

One of the shortcomings of the discrete mass representations is that it 
tends to distort the low degree harmonics. Since short arc solutions are only 
slightly sensitive to these harmonics, it was decided to impose constraints on 
the solution such that the low degree harmonics would be invariant from the 
GEM-6 solution. We set the limit at 4, i.e., the mass layer obtained from 
solving (1) should contain no harmonics of degree 4 or less nor zonal harmonics 
through degree 15. The constraint equations may be expressed as: 


362 


362 


S Am.KT (4>.\.) = T, - m. kT (4>A.) 
i=l ^ X 1 1 itl I 1 1 


1 = 0 , 1 , 2 , 3.4 ( 2 ) 

m = 0 . . .^ 
i = 5, 6, .... 15 

m = 0 
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where the surface harmonics of degree /, order m evaluated 

at the mass latitude 4>. and longitude X.., Am. are the changes in mass from 
solving equation (1) and is the initial value of the i ^ mass. For example, 
when / = m = 0 the condition is simply 


51 Am. = - 


lO 


or 


53 Am. + 53^- ~ 0 


lO 


or that the total mass of thel.ayer is zero, thus maintaining the starting value 
of GM. When /= 1, m = 0 


53^m.z. = - 53 j^- z. 

^ 11 ^ lO 1 


the Z coordinate of the center of mass remains fixed, etc. To impose the 
constraints listed above requires 36 equations, thus effectively reducing the 
rank of the linear systems to 680. Using Morrison's method, we let x = By + r 
where B is a 716 x 680 constraint matrix, x is arranged so that the 362 masses 
precede the 354 states and the first 36 masses will be dependent parameters. 
We solved the reduced system: 


B 'NBy = B'(A'Wb - Nr) 


(3) 


The constraint (2) may be expressed in matrix rotation as 


(C, C^) 


'X, 


= a 


(4) 


Cj^ is an invertible 36 x 36 square matrix, is 36 x 326, o is 36 x 1; C 2 is 
just the right side of equation (2). Solving for Xj from (4), 


Xj = C’No- - C^ x^) 


Identifying as the independent variables y, the matrix B has the form 


f-^'l'^2 A All 


B = 


V 


B, 


21 

V “ 


oA 
0 


B 


32 


J 
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where has dimensions of 36 x 326, is the identity of 

dimension 326 and B^^ i® the identity of dimension 354, r is a column vector 
whose first 36 elements equal a and the remaining elements are zeroes. 
Equation (3) can now be solved for y from which x follows. 

The method was applied to Mariner 9 analysis and was revived for the 
present application due to the similarity of the problems. 

2.9 Solutions to the Mass Layer 

-12 

The starting value of the masses were set at 10 earth masses so as 
not to disturb the already converged state vectors. Using GEM-6 as the base 
model, the RMS residuals for the various arcs ranged from . 0055 to , 013 Hz 
with an average of .0094 Hz before the mass solution. The very low residual 
level indicates that the masses should be small (nominal noise level is . 0075 Hz 
and this was the input value for data weights). While much larger anomaly 
induced signatures could be obtained by using a lower order (say 2,0 or 6, 6) 
model as base, thus resulting in larger mass values, we chose GEM-6 to 
start from for the following reasons: 

• We wish to evaluate the amount of signal not represented by 
GEM-6 

• Since so much of the anomalous gravity power is already in the 
base model, the initial states will be only very slightly corrupted 
by model errors and, hence, closer to the correct solution. 

This means that the functions should be more linear and fewer 
iterations are required in comparison to a simpler base model. 

Our experience strongly dictates that a single iteration is all 
that's necessary. 

The solution to the mass values were not uniquely determined because 
of incomplete coverage and the weak signal level. The mass magnitudes 
depended on the assigned a priori standard deviation, (j . As the a priori 
value is raised, the total RMS mass increased while the RMS data fit 
decreased. The variations are illustrated in Fig. 5. 
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Figure 5. RMS Mass and RMS Doppler Residual vs 9 

3l 


A decision had to be made on the level of a . One would like to obtain 

3 . 

a small RMS mass along with small residuals. The asymptotic behavior of both 

curves in Fig. 5 indicates that the a priori should be set in the range of about 
- 7 - 7 

. 3 X 10 to . 5 X 10 earth mass. Examination of the gravity disturbance 
plots of Figs. 6 and 7 reveals that the acceleration patterns are consistent 
i.e. , the highs, lows, and zeroes all remain in place, but the magnitudes 
change. Furthermore, over the regions of good coverage, the actual distur- 
bances vary considerably less than the ratio of .3 to . 5 which is also encourag- 
ing. Fig. 5 shows that the RMS doppler residual decreases almost exponentially. 
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“ 7 

At a = . 3 X 10 earth mass, 63,2 percent (1 - l/€) of the total decrease 

^ -7 

had been achieved. A further increase of a to . 5 x 10 em attains a net 

a 

decrease in RMS of 13 percent (of the total change) but increases the RMS 
mass by almost 45 percent with most of the large changes outside the regions 

_ 7 

of good coverage. We therefore adopted the case with = . 3 x 10 as 
our best solution and ascribe an uncertainty corresponding to the differences 
between solutions. Over the regions of good coverage, the gravity distur- 
bances differ by about . 3 to . 5 mgals. 

Some statistics on the final solution are: 

Number of arcs 59 

Number of observations 4310 

RMS doppler residual (pre) .0094 

RMS doppler residual (post) .0070 

Mean mass -. 37 E-16 earth 

(constrained to be zero) 

RMS mass .36E-7 

Formal standard deviation of each mass . 18 to .28E-7 

_7 

Since a unit (10 earth mass) causes a maximum of 2.6 mgals at GEOS altitude, 

the gravity disturbances are on the order of 1 or 2 milligals. A tabulation of 

the magnitude and locations are given in Table 4. Each point is located at 

distance 5950 km from the geocenter. The points are geographically displayed 

in Fig. 6 where the numeral locates the surface coordinate and the values are in 
-8 

units of 10 earth mass. 

A contour plot of the vertical component of gravity disturbance at 
850 km is shown in Fig. 9 representing contributions from the mass layer 
plus GEM-6. The contour interval is 4 mgal with K as the zero level. 

Since the GEM-6 contours are 10 times the contours from the mass layer, 

Fig. 9 is very close to the contours of GEM -6 by itself. 

It should be mentioned that what is plotted are accelerations and not 
gravity anomalies. A positive or mass, excess causes a negative acceleration, 
i.e. , an acceleration along the negative radial direction. (This is opposite in 
sign to the usual anomalies where a positive mass causes a positive anomaly. ) 
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Table 4. Mass Locations and Values from SST 



LAT 

LONG 

MASS 




1 

-9C .C 

C.O 

.10366-07 

.17426-07 

53 

-tl.9 

121.3 

2 

-79.4 

r,.<i 

54 

-61.9 

-166.7 

3 

-72.6 

irs.o 

-.48826-07 

55 

-61.9 

-94. 7 

4 

-68.9 

-72.0 

.14C0E-07 

56 

-61 .9 

-22. 7 

5 

-61.9 

49.3 

.43146-07 

57 

-61.9 

22.7 

6 

-58.3 

-144.0 

-.25596-07 

58 

-61.9 

94.7 

7 

-52.6 

-1C 8.0 

-.22436-07 

59 

-61,9 

166.7 

8 

-51.6 

-127.0 

-.38756-07 

5C 

-61.9 

-121.3 


-42 *u 
-40. y 
-37.1 
-31.7 
-31.1 
-29.4 
-26.6 
-21.5 
- 20.1 
-p.9 
- 10.6 

-.9 

S:3 

9.2 

10.1 


.0 - 


icr.a 

- 

1P5.0 

-120.4 
-1*5.6 
-> 2.0 

175.2 
-126.4 

150.5 

36.0 

114.2 

54.0 
120.4 
-36.0 

132.2 

-76.8 
9 


96.5 
17.4 

leoT^ 


m 

29.4 

m 

40.0 

th 

ii:! 

60.9 -36.0 

72.6 146.0 - 

79.4 36.0 


)S.9 
2.6 
► 9.4 
-79.4 
->9.4 
-79.6 
-79.6 
-72.6 
-72.6 
-72.6 
-72. 


-96.8 

18C.0 

ilhi 

-36. 
146. 
36. 

ili:i 

•144.0 
-72.0 
36.0 
180.0 
-1C8.I) 
-3< 


-72.6 -36.0 - 

:|8:l »?:8 : 


-68.9 

-60.9 


144.0 

•144.0 


.3912E-07 
.2147F-07 

.4005F-07 
.4739F-07 
•.9788F-O0 
•.1469F-07 
.3834F-07 
.6849F-O0 

.5238F-08 

.8294F-09 

.3580F-07 

.4944F-07 

•.3335E-07 

.1051F-07 

•.3600F-07 

.1324F-06 

-.6847F-08 

:imiU 

•6240F-08 
.2193F-07 
-.5209F-07 

.3727F-07 

-.1250F-07 

.6093E-07 

.1092F-07 

limKU 

•.4587E-ilT 

.1494E-07 


- 


63 

64 

65 

66 

67 

68 
69 
7C 

?l 

74 

75 

76 

77 

78 

79 


n : 


4c 

43 

44 

45 

46 

47 
38 
49 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

nc 

193 

134 


6i .9 
58.6 
58.3 
-58.3 
-38.3 
-52.6 
-32.6 
-52.6 
-52.6 
•51.6 
•51.6 
-31.6 
-51.6 
•51.6 
-fl.6 
•51.6 
•51.6 

47.7 
47.7 
47.7 
-42.0 
-42.0 
-42.0 
-42.0 
-42.0 
-42.0 
-42.3 
-42.0 
-42. f 


-49.3 

L.C 

72.0 
144. C 
-72. C 

36. C 

18J:P 

-36. C 

17.0 
89. C 

161. C 
-55. t 
53. C 
127. C 

"i61.w 

-89. C 
-17. C 
( .«. 
144. C 
-144. C 
-72. i 
43.2 
115.2 
-172.8 

-100. e 
-26.8 
26.8 
172. 


-3 


4C.0 

-37.1 

■37.1 


8 

-115.2 
-43. 2 
85.2 
.5' 
l3I 
-58.8 
58. 8 
130. 8 
-157.2 
-85. 2 
-13. 2 

i-U 

-72.0 

-144.0 


.1222F-07 
.4745F-07 
.2C76F-07 
.1375E-07 
.3267E-67 
.2e38E-C7 
.8543E-Ct 
-.8579E-0fc 
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Table 4, Mass Locations and 



from SST (Continued) 
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.2153=-:7 
' • 13 m 7 1 ~C 7 
.3136E-07 
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.3286E-07 
.i265E-j7 


139.2 

-168.9 

- 6.8 

6.8 

76.8 

168.8 

-139.2 

-67.2 

36.0 


-1C8.0 

-36.0 

119.6 

-26.5 

26.5 
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168. § 
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- 12.6 

86.6 
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.1352E-C7 
-.7660E-07 
.2220E-Q7 
.6053E-07 
-.2950E-U7 
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.2636E-J7 
.3311E-07 
-.6935E-C8 
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-.6297E-37 
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Table 4. Mass Locations and Values from SST (Continued) 
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52.6 

92.6 

tI.E 

.2955E-07 

.1C66E-07 

399 

350 
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68.9 
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36.0 
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-.25085-07 

333 

52.6 

199. C 

.9978E-08 
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68.9 

138.0 


52.6 

-199.0 

-.9967E-07 
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68.9 

16J.G 

-lOS.O 
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52.6 

-72. J 

-.57335-07 
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68.9 

-.5257E-07 

I3» 

58.3 


-.5973E-07 
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72.0 
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61.9 
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K = 0; Each level = 4 milligals, L = 4, M = 8, etc. 


Figure 9,. Total Vertical Acceleration Sensed by GE0S (GEM-6 plus mass layer) 



Also we did not Include a Brun's term. The same convention was used in 
plotting both Figs. 6 and 7. The positive contours correspond to regions of 
mass deficiency or negative surface anomaly. Some significant features 
of the acceleration map are: 

Negative Accelerations (Positive gravity and mass anomalies) 


Latitude 

Longitude 

Acceleration 

Total Mas s 
in Vicinity 

10 

20 

-1.2 

1. 6 X 10“'^ 

46 

9 

-1.6 

+2 

21 

-28 

-1. 6 

2 

Positive Acceleration (Ne 

gative gravity and 

mass anomalie 

61 

-15 

1, 2 

-2 

21 

-2. 5 

1.2 + 

-3 

Doppler 

Residuals 




Some typical residuals before and after the fit are shown in Figs. 10 
and 11, the latter improvement is less evident than the former. To show 
the effects of various levels of modeling, we also included residuals relative 
to a 2, 0 earth model, GEM-6 truncated at the 6^^ degree, full GEM- 6, and 
GEM-6 plus mass layer all for rev 268. Evidently, improvements are difficult 
to come by after GEM-6. The RMS levels were: 


Model 

RMS Doppler Residual 

Principal oblateness only (2, 0) 

. 075 Hz 

GEM-6 terms to 6-6 only 

.020 

GEM-6 complete 

. 0070 

GEM-6 plus mass layer 

. 0055 


Figure 12 shows residuals of a pass not included in the least squares solution. 

It can be seen that the irnprovement (RMS = . 008 for GEM-6 vs . 0057 for GEM-6 
plus mass layer) is on the same order as the passes illustrated above. 
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DOPPLER. H* 




1330.00 1344.00 1362.00 1300.00 1368.00 1376.00 1384.00 1392.00 


TIMC> mlnulaft 


A 


B 


C 


D 


A. ^ 4 Q RMS = . 075 h*; B. GEM-6 through degree 6 RMS = . 02 hz; 

C. Full gem- 6, RMS = . 007 hz; D. GEM-6 plus mass layer, RMS = . 0055 hz. 


Figure 10. Doppler Residuals for Rev 268 
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Figure 12. Doppler Residuals for Rev 269 



The very low level of signature to begin with coupled with the excellent 
starting states and the highly linear nature of the method enabled us to obtain 
convergence in a single iteration. This was evident in evaluating residuals 
after the solution. We therefore did not reevaluate the partials for further 
iterations. 

2. 1 1 Conclusions and Discussions 

The principal conclusions from our study of the SST data are: 

• SST tracking is a viable technique to measure the earth's near 
field. Excellent coverage can be obtained from just a few 
stations with little or no tracking gaps. 

• At GEOS-C altitude, the residual accelerations relative to GEM-6 
are on the order of 1-2 milligals. With a measurement noise of 

. 0055 - . 007 Hz observed on many of the doppler passes, the 
accuracy of the short arc method of inference is about . 3 to . 5 
milligals over the regions with reasonably good coverage. 

• The advantages of the short arc dynamical method are: high 

degree of linearity, short integration time, little aliasing due to 
plant noise, and adaptability to incomplete data coverage. The 
method is a reasonable approach for the analysis of SST data for 
the recovery of gravity disturbances at altitude. 

• Principal deficiencies in the current solution are: 

a. Degraded geometry for many of the passes, high sensitivity 
passes with small GEOS to ATS angles account for only 

/ 

75 percent of the data. Even for the ones with good 
geometry, sensitivity at the beginning and end are reduced. 

b. Much of the earth's surface is not covered by good SST 
passes. 

c. The GEOS orbit altitude is excessively high to sense some 
of the shorter wavelength terms. 

d. Many good geometry passes contained dropouts or noise 
bursts at small GEOS to ATS angles. Even the nominal 
noise level is too high to sense the residuals relative to 
GEM-6. An improvement by a factor of 2 or 5 would be 
highly desirable. 
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2 . 12 


Recommendation 


Asa project for the future, we would suggest the following experiment 
to map the remainding earth field. Launch a low altitude, polar satellite with 
perigee / apogee altitude of approximately 200/ 1000 km and track it with SST 
doppler from a high altitude polar satellite. The advantage of such a configu- 
ration is that both poles as well as the equatorial zones will be completely 
covered by passes with near ideal geometry. All of the non-hardware related 
deficiencies in the present data can be removed. All intermediate altitudes 
can be sampled as the argument of perigee precesses. 


III. Altimeter Analysis 


The objectives of the analysis are: 1) derive smoothed estimates of 

geoid height from multiple passes over a given area, 2) obtain a representa- 
tion of the geoid by a set of source parameters more compact than the 
original data set. These sources can be used to extrapolate the geoid to 
local areas not actually covered by measurements, and 3) infer gravity 
anomalies from the altimetry by virtue of the causal relations between 
geopotential, mass, and gravity. 

In view of the very extensive amounts of altimetry from GEOS, our 
scope is necessarily limited. We were designated as one' of the investigators 
to analyze data over the Atlantic Test Area. Our initial plan was to independ- 
ently determine GEOS ephemeris from C-band laser and other ground tracking 
but it became apparent very quickly that the available resources do not permit 
such a course. 

It was our intent to separate the geoid undulations into low and high 
degree components in the same manner as the SST analysis of Section II. 

The motivation here was the same as before--to see what part of the geoid 
was represented by the low degree harmonics and determine the spacing of 
a density layer for representing the remainder. Once again GEM-6 served 
as the base harmonic model. 

An outline of our tasks, then, is as follows: 

t 

(1) Select data base from altimeter tapes 

(2) Edit data to remove outliers and land mass reflections 

(3) Apply tide and known bias corrections 

(4) Make crossover adjustment on heights 

(5) Establish grid of representation elements 

(6) Calculate, geoid residuals relative to GEM-6 

(7) Infer source magnitudes from geoid residuals of step 6 

(8) Calculate gravity anomalies from inferred sources and 
compare with observations 
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Data Base 


3. 1 

The original data base consisted of 14 CALTOR BCD tapes containing 
approximately 150 passes. An initial screening chose passes contained within 
the quadrilateral shown in Figure 13. The figure is approximately equivalent 
to the Mission Plan test area in size with some displacements in the vertices. 
All points interior to the quad were accepted and the remainder rejected. 

The reason for displacing the vertices was to choose sides which were very 
nearly parallel to GEOS passes so that each one would be essentially complete 
except for possible truncation at the end thus assuring the maximum number 
of orbit crossovers. From the 140 original passes, 82 were initially selected 
with approximately 7000 observations. These had to be edited for outliers 
and land mass reflections of which there were a considerable number, because 
the southwest edge of the quad contained many islands. Table 5 lists the 
passes actually used in the analysis. The final base contained 74 passes with 
~«6050 data points at approximately 3 second intervals. 

3.2 Cros sing Adjust 

After the edit, the next task was to adjust the height of the passes for 
possible biases and satellite radius errors. The adjustment was made in a 
least squares sense by minimizing the quadratic sum of the geoid height 
differences at the crossover points, i. e, , we find corrections, ^Rj^, 3R^ which 
minimize S where ’ 



m 




(N. - 


N. + 
J 


<5R. 



m = number of crossovers 


N., N. are the geoid heights (sea height less tide) at the crossing of the 

th \h ^ 2 2 

i*" and j arc. An a priori reciprocal variance of (7 = 100 m was added to 

the diagonal of the normal equations to stabilize the solution. The adjustment 

should converge in one iteration since the partial derivatives are constants 

equal to + 1. A ctually we made two iterations because passes were deleted 

after the first correction due to some rather large residual differences even 

after the biases were applied. Some statistics on the crossover adjusts are: 
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Table 5. Summary of Altimetry Passes 
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DT 
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AO 
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55 

A 

Z 
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117 
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2.11 
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A 
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30719 
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91 

A 

H 
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C 
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A 
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211 
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199 

5.80 

90 

0 

9 
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70 

-2. 17 

32 

A 

10 

33 9 

123 
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217 

5. 89 

95 

C 

11 

374 
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35632 
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-5. 86 

30... 

A 
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843 

159 

48711 

184 
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91 

A 

13 

851 

1£3 

13219 

64 

34. 47 

29 

A ^ 
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857 

160 

47889 

194 
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90 

A 

15 

110 7 
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. £5 

91 

C 

16 

1156 

161 

5 9347 

127 

1.81 

36 

A 

17 

1505 

2C6 

32719 

175 

56 

67 

C 

18 

1511 

206 

67514 

186 

3.42 

92 

A 

19 

158 2 


69146 

182 

7.72 

9C 

A 

70 

1590 

212 

33515 

219 

-9. 99 

ICC 

Z* 

71 

1647 

216 

35975 

209 

-8.46 

94 

c 

22 
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74053 
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1.91 

93 

A 

23 

1875 

223 
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90 

C* 

2<* 
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24 3 
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-3. 93 

90 

c 

25 

2037 

2-,3 
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-1. 45 

93 

A 

26 
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76509 

213 

5 . 38 

99 

A 

27 
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85499 
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89 

A 

28 

2102 
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49679 
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77 

A 

29 

21C8 

2 48 
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1.11 

90 

A 

30 

2145 

252 

52343 

186 

-14. 88 

92 

Q •» 

31 

2159 

253 

515G7 

207 

.33 

96 

c 

32 

2244 

258 

53135 

213 

2. U 

101 

C 

33 

530 

137 

38C76 

180 

-3. 62 

65 

A 

3<i 

587 

141 
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166 

-8.32 

82 

A 


35 

2372 

267 

5 7224 

199 
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98 

C 

36 

235 0 

266 
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205 
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A 

37 

2358 
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58061 
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89 

C 

38 
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34837 
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-4 . 5 C 

84 

A 

tC 
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113 
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1.16 

59 

C 

40 
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71 

c 
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5C 
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A 
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57 

A 

45 
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2.45 

88 

A 

46 

1960 
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93 

C 

47 

1965 
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81411 

176 

-. 57 

85 

A 

48 

1974 
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45 78 9 

198 
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9 2 

0 

49 
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<<.2 

49091 
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79 

c 

5C 


2s7 

50 715 

187 

•4 . C C 

92 

c 

51 

2770 

29p 

6 8661 
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3. 11 
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c 

52 

2890 

3C4 

21576 

186 

1. 13 

91 

A 

53 

2898 

3C4 
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1C2 

C 

54 

2933 
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2525G 
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69 

t 

55 

2941 
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c 

56 

3C69 
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30 142 
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65 

c 

57 

3097 
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-. £8 

ICC 

c 

58 

3154 
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c 

r c 
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A 
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45 

A 

61 
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3.1c 

86 

A 

62 

2699 
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98 

Z 

c3 

2756 
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6 54 9 6 

198 

3.26 

98 

c 

£4 

215 0 
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15 = 1 
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. £ 0 

94 
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65 

230 7 

<£3 

3996 
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-3. 15 

89 
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66 

53 8 
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2443 

199 

-2. 29 

95 

C 

67 
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159 

1393 0 

182 

-.CO 

9G 

c 

68 
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22936 
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97 

c 

69 
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186 

60974 
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- . 76 

97 

A 

70 

156 8 
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6 59o8 

181 

-1. £1 

86 

A 

71 

595 

ls2 

49j5 

2G9 

-3.16 

95 

C 

7 2 

1576 

211 

3^3^7 
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-1.56 

94 

0 

73 

452576154 

70 974 
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• 41 

58 

C 

74 

458276158 

734 30 
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-1. £4 

£0 

A 


• 
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NP = pass number; All passes in 1975 except 73, 74 which are in 1976; DT = duration of pass; 
BIAS = best estimate of bias from crossing adjust; NPT = number of observations; AD, A = 
ascending, D = descending. 
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Case 2 



Number of Passes 82 

Number of Intersects 1448 

RMS Height Difference - meters 

Before 1 6. 0 

After 1.37 

RMS Bias - meters 11,26 

Mean Bias - meters . 031 


74 

1148 

1. 28 
.92 

6.04 (3.3) 

. 28 


The corrections obtained from the second case were added to the observed 
heights to obtain the corrected geoid. If five of the final 74 biases with 
magnitudes greater than 10 meters were excluded, the RMS bias was reduced 
to about 3, 3 meters. The probable accuracy of the ephemeris radius is then 
3. 3/v/2 Si 2.3 meters. One bias which cannot be removed by the preceding 
is the absolute or constant calibration bias [11]. That one was taken out just 
prior to the mass layer regression. 

Before proceeding to the next steps of the analysis, we make a digression 
on the problem of representation. 


3, 3 Representation 


Some candidate elements of representation are: 


(N = geoid height, g = gravity anomaly, m = mass element) 



(1) 

N —Ag 

Molodenskii 1962 fl23 


(2) 

N ^m 

Chovitz [ 13 ] 


(3) 

Ag-* N 

Stokes' formula Heiskanen, Moritz [l4] 


(4) 

m N 

Classical potential 


(5) 

N -f(<|>,\) 

Polynomial expansions in latitude (4>) and longitude (\) 
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The first two express the sources as a convolution of the measured 
height and an influence or kernel function over a sphere. They require obser- 
vation data over the entire surface. Since this is generally not available, the 
currently popular approach is to use observations near the point of interest 
and substitute geoid values from a spherical harmonic expansion outside of 
an assigned polar angle limit. The same requirement exists for methods 3 
and 4. Our main interest is centered on method 4. The others are being 
pursued by other investigators. 

There is considerable similarity between Stokes' formula and the classical 
potential 

^ ~ 4iTYa / S(4<) Ag(4i) d(7 Stokes 

(7 


N = a^ /* = a^/* Classical Potential 

where S(4') is Stokes' kernel, Y and a are mean gravity and earth radius, 
respectively; d<7 is an element of area on the surface c, P is the distance from 
the source to the field point; X(4') is an equivalent surface density, dm is an 
element of mass (in units of earth mass). The representation elements are 
Ag, da and Xda, respectively. Hence, the difference in the two equations is. a 
matter of the difference between Stokes' and the 1/P kernel. Assuming that 
the density elements are located on the surface of a sphere and letting 

t - 

2a _ 
p ~ t 

S(‘|^) = ~ - 6t + 1 - 5 COS 4* - 3 cos 4 ln(t + t^) 


For small values of both kernels have a 1/t type singularity; has an extra 

• 1 ' ^ C\ 

factor of 2 in the denominator . Figure 14 shows the two functions for 4> ^ 10® 


I 


5*C 

The extra factor of 1/2 is immaterial since it cancels out in the ratio P(4')/ 
p(4* < 1°). This ratio is indicative of how rapidly the integral converges with 4. 
Hence, it is simpler to compare Stokes' with esc 4/2 rather than (1/2) esc 4/2. 
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Figure 14. Stokes' and Inverse Distance Kernel 
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Magnitudewise, Stokes' kernel is larger over most of the range but it passes 
through, zero at two points with two sign changes. The fact that the 1/P kernel 
has the same sign is of little consequence since the source elements are 
certain to change signs many times. On the other hand, the rapidity of the 
magnitude decay dictates the rate of convergence with increasing 4'. For 
practical purposes, the two representations are essentially equivalent since 
the kernels behave similarly near the origin and converge at approximately 
equal rates. 

Equating the two expressions for N and retaining only the critical region 
near 4* = 0, one obtains: 


1_ 

4TrYa 


SAg da = 


Z da 


or 


- 4 ^ = 

4TrYt 


"2T 


which is the result of applying Gauss' theorem V • (Ag) = 4 itGMX to a large 
plate. 

Since both kernels become singular as 4'— *0 even though the integrals 
are well behaved, some care is needed to average and weight the elements 
with small 4. One simple way to circumvent the singularity is to remove the 
discrete elements from the surface. Another is the use of finite (or distributed) 
m.asses either at the surface or at depth. Both methods result in a smoother 
variation of the potential (or geoid) than surface sources and are likely to be 
better representations of the actual data. 

f 

While the finite distributions are intuitively more appealing than point 
sources, the computation of the potential, even for very simple shapes is much 
more complex. To get some feel for the differences between point and dis- 
tributed masses, we calculated the potential of a buried disk at various depths 
and compared its potential with that of a corresponding point mass. The 
potential of the disk was obtained by numerical quadrature (see Appendix I 
for algorithm). Two cases were examined. 

a. Disk and point mass at the same depth but with different masses 
such that both induced a geoid perturbation of 1 meter at a point 
on the surface directly over the mass (i. e. , on the axis of the 
disk). 
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b. Total disk mass equals point mass but with the latter at a greater 
depth such that both induce perturbations of one meter at the 
surface point directly over the mass. 

We considered depths of b and 2b where b is the radius of the disk. For the 
sake of brevity, only the results of (a) will be shown since differences for (b) 
are smaller. 

Figure 15a shows the geoid heights N^, for the disk and point source, 
respectively, versus normalized distance from the axis P = r/b, and ^ = Z/b 
with b - 20 nmi where Z and r are coordinates along and perpendicular to the 
disk axis, respectively. The percent differences between the two are shown 
in Figure 15b. Significant differences exist mainly for small values of P as 
one expects. When these normalized variables exceed 2, the disk begins to 
resemble a point source. In the critical region of small P, ^ , the differences 
grow almost linearly from zero to the steady state values of ~18 and 5 percent, 
respectively, with the disk having a slower decay rate. 

The gravity anomaly and their differences shown in Figure 16 are another 
matter. The fall off of gravity for the point mass is significantly sharper 
than the disk. The combination of both observables, i. e. , gravity and geoid 
(or potential) could eventually provide a basis for the choice of a model. Since 
the current altimetry data base completely overwhelms the amount of available 
gravity data, we decided to start with a point mass model because of its 
simplicity and rather mild geoid differences noted above. If and when the 
residual signatures in altimetry or gravity dictate a more refined model, 
software for computing the potential of the disk or even solids will be available. 

Having chosen the point mass model, there is still the matter of a starting 
grid spacing and depth. Again, there is no hard or fast rule. We chose an 
initial spacing of 111 km (1 deg) and depths of 55. 5 kna and 111, 1 km. Experi- 
ments will be the guide in dictating other choices. 

3. 4 Some Remarks on Representation 

The basic problems in analyzing GEOS altimetry are representation and 
extrapolation to unmeasured points. The inference of gravity anomalies and 
disturbances is a related third problem. Given adequate coverage, the geoid 
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Figure 15 , Geold Difference Between Point and Disk Mass 
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Figure 16, Gravity Disturbance Difference Between Point and Disk Mass 
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could be approximated locally by polynomials which can be spHned to obtain 
any degree of continuity fis]. To be practical, the expansions must converge 
rapidly. The fitting functions are expected to be valid extrapolations of the 
geoid at unsampled points in the vicinity of the measurements because the 
potential is smoothly varying. (At least over distances on the order of a 
correlation length.) We also expect the derivatives of the fitting functions to 
be reasonable estimates of the vertical deflections and local gravity disturbances 
with errors which depend on the extent of the measured region. 

Use of the fitting functions to extrapolate potential or gradient to another 
altitude is related to the third problem. Classical integral formulas allow 
such an extrapolation given coverage over a closed surface. In treating local 
areas, one justifies the extrapolation by noting that the kernels of the integrals 
fall off rapidly with distance and one can truncate the integration at some 
"small" distance within acceptable error limits. Physically, it is a matter of 
causalty, local mass anomalies induce local disturbances. 

Some desirable properties of the point mass kernel at depth are: 

(1) Function is harmonic and simple to compute 

(2) Can be fitted to match boundary conditions at surface 

(3) It is smoother than surface mass kernel with surface sources 

(4) Has the correct boundary condition at infinity 

(5) Has some physical significance 

(6) May be used to represent or infer gravity 

3. 5 Truncated Kernels 

One aspect of the gravity anomaly to geoid (Ag — ♦ N) problem which has 
received some attention in recent years is that of missing measurements 
beyond some limiting angle away from the field point. In addition to the 

substitution of anomalies from harmonic models in the region iIj > i{j , 
suggestions have been made to uSe a so-called "truncated" kernel to further 
lessen the error caused by incomplete data coverage. Stokes' kernel has been 
analyzed at some length in this respect. The idea is to remove the low degree 
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terrns in the Legendre series expansion of the kernel and use orthogonality 
properties to show the improvemeiits thus derived. 


The same type of removal could be done in the case of the point mass 
kernel, i.e. , in place of a/P, we used the kernel where 

ji , . n+ 1 


P (cos 4*) 
n 


Rather than go through the lengthy derivation of the error estimates, we 
graphically illustrate the functions K^. The behavior of the for I ~ Q, 4, 8, IE 
and 2 < 4^ < 65 deg”'' are shown in Figure 17. The more rapidly Kj^ decreases 
the faster the convergence with 4^, Using the point 4^ k . 27^ as reference where 
f(. 27°) “ loo, the values of for 4* > 2. 5 deg are all less than 15. In the 
region 2. 5 < 4^ < 30 degrees, the modified kernels are considerably smaller 
than the one with / = 0. Since all the through / 12 are within 12 percent 

of Kq for 4 j < . 27°, the more rapid decrease from about 2.5 to 30 degrees 
implies less error arising from loss of data within this zone. 

3. 6 Mass Layer Regression 

A uniform grid with spacing of one degree between points was located 
beneath the quadrilateral in Figure 13. Inference of mass magnitudes at these 
points follows the usual weighted least squares algorithm. A single iteration 
suffices in each case because of the small residual amplitude relative to GEM“6. 
No a priori was added to the diagonals because of good data quality but the 

t 

boundary points might have benefited froin some degree of stabilization. Two 
mass depths were considered -- -55. 5 km and -111. 1 km. 

The coinputed geoid height from GEM- 6 followed Rapp's algorithm [l 6] 
with flattening f s 1/298. 255. The term was essentially zero while 

was altered from the GEM-6 value of . 536 x 10 ^ to -. 254 x 10 ^ resulting in 
a geoid shift of about -4. 25 meters over the quad. 


>}< 


a (6378, - 55. 5) km, r = 6378 km, i. e. , masses at 55.5 km depth. 
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Figure 17. Modified Point Mass Geoid Kernel vs ^ 
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Table 6.. Mass Solution from Altimetry, Depth = -55, 5 km 
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Results 


3. 7 

The smoothed geoid contours due to the mass layer over the test region 

_9 

are shown in Figure 18 and the inferred mass distribution in units of 10 
earth mass at depth -55. 5 km are shown in Figure 19. Some statistics for 
the best solution (listed in Table VI) and several related ones are given in 
Table VII. 


Table 7. Statistics of Mass Layer Solutions 



<0 = «• 

536 X 10"^ 

^40 = 

-0. 254 X 10"^ 


No Bias 

With Bias ' 

No Bias 

With Bias*’' 

Mean residuals GEM-6 only (m) 

1. 556 

5. 106 

-2. 691 

. 858 

RMS residuals GEM- 6 only (m) 

2. 967 

5. 697 

4. 293 

3. 453 

Mean residual (with masses) (m) 


-.017 

-.018 

-. 017 

RMS residual (with masses) (m) 


. 717 

. 707 

.711 

-9 

Mean mass (10 earth mass) 


. 544 

-. 116 

. 196 

“ 9 

RMS mass (10 earth mass) 


4. 458 

4. 218 

4.330 


The smallest mean geoid residual (without the mass layer) is obtained for the 
case = -0. 254 X 10 . Calibration-bias-corrected residuals for this case, 

with and without the mass layer, are shown in Figures 20a and 20b, respectively 
Many of the preadjust residuals are in the 6- to 7*- meter range. Figure 20b 
indicates that the noise level of the data is even better than the RMS of 0, 71 
meters because there are still clearly discernable trends. Some of these may 
be due to ephemeris errors. It appears that the one-degree spaced grid is 
adequately dense for most of the regions. We have not attempted to insert more 
masses, even though some densification is evidently called for in some areas. 

Since each observation was weighted by the inverse standard deviation, 

the mass estimates are also minimum variance. The formal standard deviation 

-9 -9 

of each mass ranged from about 3 x 10 earth mass at the edges to 0* 3 x 10 

earth mass near the center. Truncation effects have not been included and 
would increase the uncertainty much more, particularly near the boundaries. 


‘Calibration bias = -3. 55 meters. 
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Fig 8A, K = 0, each level equals 1 meter, L = 1, M = 2, etc. 
Fig 8B, U = 0, each level equals 3 meters, T = -3, S - -6, etc. 


Figure 18. Geoid Heights from Mass Layer, Total Geoid 
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Solutions were also obtained for a mass layer at a depth of -111 km 
and for cases using a kernel with low degree harmonics removed. The primary 
effect of doubling the depth is an increase in RMS mass by a factor of 4 as one 
might expect. The much larger fluctuations indicate that the correlation 
between neighboring masses becomes excessive at -111 km. The fit of resid- 
uals decreased only about 2 percent when the depth was doubled. A truncated 
kernel caused significant mass changes mainly near the edges and hence a shift 
in the mean. RMS residuals remain essentially unaltered. 

The contours of Figure 18a are mostly gentle. A negative region with 
trough of -5 m exists at about 33N, 285E. Highs of 10 meters can be seen at 
+27. 5*^N, 292°E and along the southern border. Some of these are probably 
due to edge effects. The total geoid in Figure 18b is the sum of GEM-6 and the 
mass layer contribution from Figure 18a. The composite model represents 
the smoothed observations to an RMS level of , 71 m. 

3. 8 Gravity Anomaly Differences 

The Defense Mapping Agency provided a total of 67 1° x 1° anomalies 
in the region 18° to 40° north latitude and 27 5° to 295° longitude. Unfortu- 
nately, only 11 of them were interior to our quad and, of these, most were 
near the edges where the accuracy is degraded (see Figure 9)* Also, the 
average standard error of the observations was 13. 6 mgal. The mean and 
RMS difference of the 11 interior anomalies with respect to the combined 

model was approximately +20 mgal independent of whether the mass layer 

* 

was at a depth of -55. 5 or -111 km' . Even though there were some rather 
large residual changes for individual points, the mean and RMS did not depend 
on whether the computed values were point anomalies or anomalies averaged 
over the corners of the square. Nor do they depend on whether the full or 
truncated kernel propagated the mass layer. The mean difference decreased 
2. 5 mgal to 17. 5 mgal when the calibration bias was removed ' and the RMS 
decrease was somewhat smaller at 2.2 mgal. The remaining bias discrepancy 
cannot be readily explained at this time. The gravity sample is really too 
small and the observation uncertainties too large. It should be mentioned that 
in all of the cases, the RMS difference about the mean was only 7 milligals. 


These results can be shown analytically. 
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Table 8, Gravity Anomaly Differences 


+ 


A 

g 

Se 

/V> 

8h 

8h 

dg 

8h 

dg 

34. 5 

285. 5 

-49. 0 

-27. 3 

-21. 7 

-26. 3 

4. 6 

-41. 6 

19.9 

32. 5 

281. 5 

- 5.0 

-19.0 

14. 0 

-8. 8 

22. 8 

7. 5 

6. 5 

31. 5 

283. 5 

-40. 0 

-25. 9 

-14. 1 

-42. 8 

28. 7 

-25. 7 

11* 6 

31. 5 

290. 5' 

-35. 0 

-44. 8 

9. 8 

3. 5 

6.4 

-4.8 

14. 6 

31. 5 

291. 5 

-26., 0 

-46. 3 

20. 3 

-2.4 

22. 6 

-2. 3 

22. 5 

30. 5 

291. 5 

-38. 0 

-48. 5 

10. 5 

-.5 

11. 0 

1. 8 

8. 6 

29 . 5 

290.5 

-38. 0 

-48. 7 

10. 7 

-9. 6 

20. 3 

-3. 1 

13. 7 

29. 5 

291. 5 

-31. 0 

-50. 3 

19. 3 

19.0 

. 3 

3.3 

15.9 

27. 5 

289. 5 

-30. 0 

-48. 9 

18. 9 

-5. 7 

24. 6 

-9. 6 

28. 5 

25. 5 

287. 5 

-17. 0 

-44. 4 

27. 4 

4. 1 

23.3 

8. 8 

18. 6 

25. 5 

289. 5 

-37. 0 

-49. 9 

12. 9 

-11. 1 

24. 0 

-3. 7 

16. 6 


Table 8 shows the results from an example run with the mass layer 
s!< _ 6 

at -55. 5 km, ~ -.254 x 10 and bias removed. The definitions of the 
symbols are: 

4^,^ = Latitude, longitude, respectively (degrees) 

g“ = Observed 1*^ x 1° anomaly from DMA - milligals 
g^ = GEM- 6 anomaly 

^ - gg = "observed” high degree anomaly 
gj^ = Computed gravity from mass layer 

dg = th “ ®h 

g^ - Computed gravity averaged over center and corners of 1° x 1° square 
dg = tf, - 


The large differences between g^ and g^ suggest that the gravity field due to 
the layer is too rough and that the finite shape masses should do better. We 
could weight the gravity data in the mass layer regression but the sample size 
seems too small at present to have much difference. Further investigations 
are needed to resolve the differences. Unfortunately, time and funding did not 
permit continuation of the current effort. 
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3.9 


Conclusions and Remarks 


o The crossing differences provide an excellent measure of 

consistency on the altimetry and ephemeris. With moderate 
care in editing, RMS crossing differences (after adjustment) 
of less than 0,9 meters are attainable. The RMS bias 
corrections indicate an ephemeris precision on the order of 
2. 3 meters. 

o Relative to an ellipsoid of flattening = 1/298. 255, the mean and 
RMS geoid residual with respect to GEM-6 over the Atlantic 
area examined were 0. 86 m and 3. 45 meters, respectively 
(after application of the calibration bias). The mean residual 
corresponding to the actual value of in GEM-6 is 5, 1 meters. 

o A grid of mass points with 1 degree spacing at 55. 5 km beneath 
the surface fits the residual geoid signature to a level of 0. 71 
meters RMS. Using minimum variance weighting, the formal 
error estimate for a mass near the center of the distribution 

-9 

is about 0. 3 X 10 earth mass. 

0 Thus far, comparisons between observed and inferred gravity 
have been inconclusive due to small sampling, marginal surface 
data quality, and apparent roughness of the mass layer. Further 
investigations are needed to resolve the bias of some 17. 5 

1 ..... 

milligals for the 11 points over the area. 
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APPENDIX 1 


COMPARISON OF DISK AND POINT MASS POTENTIALS AT DEPTH 


Two cashes are of interest. The first is one where both disk and point 
are placed at the same depth with masses such that they produce the same 
potential or geoid height shift at the intersection of the disk axis and the earth's 
surface. Let the depth of both masses be Z = nb where b is the disk radius 
and Z is distance along the disk axis from the center. The mass of the point 
source required to produce a height N is given by: 


m _ nb 

M " 2 ^po 

a ^ 


( 1 ) 


M is the earth's mass and a the mean earth radius. The height shift along the 
disk axis at distance ^ = Z/b | N 


'do 


N 


do 


ZTTorba 

M 


(/Th-w 


(2) 


where a- is the density (mass/area). Hence, 


cr 

M 


N 


do 


2trba 


(3) 


At points r, Z where r is distance from the disk axis (G = gravitational constant) 


N 


Gm 




2 

a 

GM 


nN 


_QO. 




(4) 


N 


I(P,C)N 


do 


2ti( 


(5) 


1*1 


where P = r/b, Z/b and I(P,C) Is an integral given in Appendix 2. Upon 
setting =1, ^ = n, we obtain the geoid shifts of the two masses at 

the same depth. Fig. 4 shows the heights for a disk and the percent differences 
defined by 


p = 100 


N , - N 
d 1 


for values of ^ = 1,2, 

The second case of interest is one in which both disk and point source have 

the same mass but are placed at different depths such that both produce the same 

2 

shift at the surface. From (1) and (4) with N =1 meter, r = 0, and m = irb , we 

V / pQ 

obtain 


Np = (6) 

Substitute cr/M from (3) and set = 1» then solve for Z = Z^, define 

d 

pa 


d = 
P 




The value of d^ from (7) is used in (4) for points off of the axis. A point mass 
at depth d^ causes the same geoid shift at the surface as a disk with equal total 
mass oriented with axis parallel to tbe surface normal and disk center at depth 
4. The percent differences for case 2 are much smaller than the previous ones 
indicating that the deeper mass approximates the disk potential better. 

Consider next the gravity disturbance, 6^, arising from the masses. At 
the same depth with masses given by (1) and (3), 


b^{P^+C^) 


(P^+t^) 
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where 6g is in units of Y's, and N is set equal to one meter. The disk 
P po . 

acceleration is shown in Appendix 2 to be (T also defined in Appendix 2): 


5gd I' 

S^d = -j=- 

2irb(yi+C -C) 

This comparison of equal masses at different depths is obtained as before by 
holding 6g^ fixed and calculating 6gp from a point mass located at ^ The 

expression is simply; 

m a^ 1 ^d 

= M b2,p2^42, = b(p2,^2, 

Note that in neither case does 6gp - 6g^ at the surface. Since the masses and 
depths were chosen to give equal geoid height shifts, they cannot be made to 
have equal 6g*s (i. e. , derivatives) also. 
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APPENDIX 2 


POTENTIAL AND ATTRACTION FROM A DISK 


The potential of a thin disk having uniform density, o', and radius, b, 
is given by (see Fig. 2-1 for the coordinates): 


V = G<r| f - - ^1- - ; S = [Z^ + - 2rr cos 6] 


f 

o o 


Using normalized variables p= r/b, C= Z/b, and integrating with respect to 
F, with the definitions: 


g = P cos 0 


,2 ^ r2,l/2 


fj = (1 - 2g + p‘' + n 
f^ = (P^ + ^2)1/2 


One finds. 


V 

Go-b 


= I 


'1 - g + £, 


1 + / .^1 “ ^2 ® ^2 • ® 


de 


( 2 . 1 ) 


V = V(P, S) were obtained by setting C at constant depth and numerically 
integrating (2. 1). 

To find the attracUoi^* 6g, we note that: 


1 

Go- 


1 iv _ _ 

Go-b a; ~ 2 " 




1 + f 


g f < (1 - g 


1 ■ 

+ f,)J 


de (2.2) 
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When P = 0, both are simply integrable, and reduces to: 



which are the elementary expressions for the potential and attraction of a disk 
along the axis. 



Figure 2-1, Disk Geometry 
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Appendix III. Error Analysis of Truncated Kernels 


Let the total potential at each field point, T, be represented by a surface 
integral over 


T - - /* jT dnri 

•{7P J I P ^ Jz ^ 


(3. 1) 


where region 1 is the part of <7 with measurements and region 2 is the part 
without. Decomposing into low and high degree components designated by 
subscripts /. and h, respectively, and using K for 1/P, the actual high degree 
part of T is given by 


Tu = T 

h 


■ " Xz " fl+2 ^ ‘*”'h 


The last equality invokes the orthogonality of and Consider the error 

incurred from using the kernels K and over the region 1 integration. For 


the kernel K, , the error is 
h 




K. dm. = /. K, dm, 

h h h h 


(3.3) 


Had one used K = + K,^, the error would be 


h’ : 


/ 

-'1+2 


K. dm. 
h h 






K, dm, 
h h 




£ 


Kl dm^ + ^ K^ ^”^h ~ S ^ 


h h 


(3. 4) 


Use of the high degree kernel reduces the error by the amount J' Yi^ 
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From equations 3. 3 and 3. 4, the variance due to the truncation are 


|a = . 


Ek = kE, K' 


respectively, 


and since |k| we can expect a smaller error from use of the high degree 

kernel. In practice, Kj^ is simply given by 

l 

= f - ^ 


n=o 


where i is the degree of removal. 
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